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In the three-dimensional Heisenberg spin glass in a random field we study the properties of the 
inherent structures that are obtained by an instantaneous cooling from infinite temperature. For 
not too large field the density of states g{iv) develops localized soft plastic modes and reaches zero 
as (for large fiels a gap appears). When we perturb the system adding a force along the softest 
mode one reaches very similar minima of the energy, separated by small barriers, that appear to be 
good candidates for classical two-level systems. 
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Supercooled liquids and amorphous solids exhibit an 
excess of low-energy excitations, compared with their 
crystalline counterparts [1], in which at low frequencies 
the density of states (DOS) g{ijj) has a Debye behavior 
gioj) (X in d spatial dimensions. This excess of 

low-frequency modes is called Boson peak ma and it is 
located at a small, but non zero frequency. 

What happens at much lower frequency? Obviously we 
find phonons, but what is there beyond phonons? Were it 
possible to disregard Goldstone bosons, a scaling g{u}) oc 
(jj^, with ^ = 3 or 4 has been suggested for disordered 
systems mis]. Still, this has not been demonstrated nor 
observed. It has been stressed by EHH] that there are 
localized plastic modes, whose spectral density reaches 
zero when uj goes to zero. These modes are subdominant 
in the small frequency region: They are called “plastic" 
because they dominate the plastic response. These extra 
small frequency modes may be related to the behaviour 
of hard spheres at jamming [US]. 

Replica theory offers an explanation for these extra 
modes. At low enough temperatures, strongly disordered 
mean field models undergo spontaneous full replica sym¬ 
metry breaking (RSB). Full RSB implies a complex en¬ 
ergy landscape with a hierarchical structure of states and 
a large amount of degenerate minima separated by small 
free energy barriers m IS]- As a consequence, zero- 
temperature equilibrium configurations can be deformed 


at essentially no energy cost through easy-deformation 
patterns, that we name soft modes [I3IS]- These modes 
are localized in space, but non-exponentially. In fact, 
the zero temperature phase transition from the replica 
symmetric phase to the RSB phase is accompanied by a 
divergence of the localization length [13- 

More often than not, finding low-lying energy minima 
of glassy systems is a NP problem [S] ■ Here we study the 
behavior of inherent structures (IS), local minima of the 
energy obtained by relaxing the system from high tem¬ 
perature (the thermal protocol should not change dras¬ 
tically the DOS, at least if we remain in the replica sym¬ 
metric phase, see the appendix). In our study we need a 
model with continuous degrees of freedom. The Heisen¬ 
berg model, where the spins are three-dimensional uni¬ 
tary vectors, is an epitome of the spin glass m- 

The global rotational-symmetry of the Heisenberg spin 
glass has far-reaching consequences. The correspond¬ 
ing symmetry in structural glasses is translation sym¬ 
metry (which has similar implications). The Goldstone 
mechanism induces soft excitations in the form of spin 
waves [20) . Even in disordered systems, spin-waves are 
efficiently labelled by their wavevector, especially at low 
frequencies (see e.g. [21]). As a consequence, we have a 
Debye spectrum g{uj) oc i.e. extended spin waves 

(for structural glasses the situation is slightly more com- 
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plicated. These symmetry-induced modes mask the 
physics we aim to investigate. 

Thus, we add a random magnetic field (RF) to wipe 
out the symmetries. Indeed, mean field suggests that, if 
small enough, our RF does not destroy the glass phase 
[26j . In a RF, spin waves have a positive frequency, even 
for vanishing wavevector (e.g. a ferromagnet in a RF has 
no soft-modes). Therefore, the RF exposes the (possibly 
localized) plastic modes that interest us. The resulting 
spectrum has no reason to be Debye as it does not result 
from plain waves. Yet, a crossover to the Debye regime 
should appear when the RF is small. A similar proce¬ 
dure of symmetry removal has been carried through in 
glass-forming liquids, by pinning a certain fraction of par¬ 
ticles UZHSI]- The Heisenberg spin offers the advantage 
of allowing us to simulate unprecedentedly large systems, 
letting us observe scalings of several orders of magnitude. 

Here, we study the ISs starting from initial random 
configurations and we do find that they are marginally 
stable states: the distribution of eigenvalues of the Hes¬ 
sian matrix stretches down to zero as a power law, and it 
is unrelated to symmetries in the system. Furthermore, 
we find that the soft modes are localized. We also take in 
account the anharmonic effects due to the complexity of 
the energy landscape. We find that the energy barriers 
along the softest mode are extremely small and that they 
connect very similar states with an strong relationship, 
that we propose as a operational definition of classical 
two-level systems (TLS) |32] . 

Model We study the three-dimensional Heisenberg 
spin glass in a RF. The dynamic variables are three- 
dimensional spins Sx, placed at the vertices a; of a cubic 
lattice of linear size L with unitary spacings. We have 
therefore N = spins, and 2N degrees of freedom (dof) 
due to the normalization constraint sj = 1. The Hamil¬ 
tonian is 

N 

^Rf{|^) — ^ ^ Jxy^x ‘ Sy ^ ^ ^x ' ^x 5 (1) 

\x-y\ = l X 

where the fields hx are random vectors chosen uniformly 
from the sphere of radius iFamp, and |s) indicates the 
full configuration of spins Sx- The RF breaks all rota¬ 
tional and translational symmetry, removing the Gold- 
stone bosons. The couplings are fixed, Gaussian dis¬ 
tributed, with Jxy = 0 and J^y = J^, where (...) is the 
average over the disorder. 

We simulated on systems of linear lattice size L = 
12,24,48,96,192. We chose always J = 1, and we com¬ 
pared it with TFamp = 0.01,0.05,0.1,0.5,1,5,10,50. In 


^ In the case of structural glasses, the authors of I22H25I find 
g{ui) ~ is valid for any spatial dimension. The origin of this 
different behaviour should be investigated carefully. 
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Figure 1. Cumulative F{X) of the spectrum of TVd. In each 
plot we show black a reference curve representing the power 
law A^'®, and a grey line indicating the Debye behavior A^'®. 
Inset'. The DOS g{oj) calculated with the method of the mo¬ 
ments. In the limit of a diagonal hamiltonian (J = 0) the 
DOS would be a delta function centered on — Hamp. This 
value is represented with vertical lines. 

the appendix we summarize the simulation parameters. 
The case iFamp = 0 will be treated in a future work [55] . 
because it requires a different type of analysis, since the 
spin waves do not hybridize with the bulk of the spec¬ 
trum. 

Density of states We calculate the dynamical matrix 
as the Hessian matrix M of Hamiltonian 0 , calculated 
at the ISs. In the appendix we report how the ISs were 
obtained, we motivate the choice of the algorithm and 
the temperature of the starting configuration for the 
energy minimizations, and we show how the Hessian ma¬ 
trix M was calculated. 

Once A4 is known, from each simulated Hamp we calcu¬ 
late the spectrum of the eigenvalues p(A) or equivalently, 
in analogy with plane waves [84] . the DOS g{uj), by defin¬ 
ing X = We measure the DOS both with the method 
of the moments |55H37], and by explicitly computing with 
Arpack |55| the lowest eigenvalues. 

We find that, although for large fields there is a gap 
in the DOS when the field is small enough the gap dis¬ 
appears and the DOS goes to zero developing soft modes 
(Fig-iii inset). 

We focus on the p{X) for small A, or even better on its 
cumulative function F{X) = p(X')dX'. If F(X) reaches 

zero as a power law, we can define three exponents S, a 


^ The method of the moments returns the full density of states, 
but it is not precise at the tails. With Arpack we can calcu¬ 
late exactly the smallest eigenvalues, but only a small number of 
them. So, when we want to show the whole spectrum we need to 
use the method of the moments, while when we show the softest 
part of the spectrum we need Arpack. 
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Figure 2. Inverse participation ratio as a function of the nor¬ 
malized rank i/2N of the eigenvector (i = 1 has the smallest 
eigenvalue, i — 2 the second smallest, ...), for -ffamp = 1 
(top) and Ha,rap = 0.1 (bottom). Inset: Correlation func¬ 
tion C(r) extracted from the eigenvectors, for fields (from top 
to bottom) Hamp = 0.01,0.05,0.1,0.5,1,5,10,50 in L = 192 
lattices. See the appendices for a close-up. For the smallest 
Hamp our data do not display an exponential decay, which 
could be caused by a localization length larger than the sys¬ 
tem size. 

and 7 , that describe how the functions g, p and F go to 
zero for small A (or w): 

p(A)^A“, F(A)~A^ (2) 

where the exponents are related by d = 2 a -|- 1 = 27 — 1 . 
In the absence of a field one expects a Debye-like behavior 
<5 = d-l = 2, a = 0.5, 7 = 1.5 HJ. 

In Fig. we show the function F{\) for fields i?amp = 
0.1,1. The plots are compared with the Debye behavior 
A^ ® and with the power law behavior A^’®, because our 
data suggests a universal behavior around 7 = 2.5 {8 = 4, 
a = 1.5) for all i/amp that does not exhibit a gap. See 
the appendix for the data on other Ffamp- 

When the field is small we remark a change of trend 
from 7 « 2.5 to 7 < 1.5 at a value A*. Very roughly 
speaking, the crossover point goes as A* ^ maybe 

indicating the presence of a boson peak. 

Localization Similarly as it happens in other types of 
disordered systems the soft modes are localized, 

meaning that the eigenvectors | 7 f„) are dominated by few 
components. A nice localization probe is the inverse par- 
ticipation ratio F„ = ■ If the eigenvector | 7 r„) 

is fully localized in one site, then V„ = 1 , whereas if it is 
fully delocalized, V„ = 1/fV. In Fig. we show that the 


® The value 7 = 2.5 is also hypothesized in [4], through a fourth- 
order expansion of a single coordinate potential of the minimum 
of the energy. 


softer the eigenvectors the more localized they are, and 
for infinitely large systems there is probably a localiza¬ 
tion threshold that separates a small fixed percentage of 
localized eigenvectors from the delocalized bulk ones. 

The localization length increases as iFamp decreases, 
see Fig. [^inset. In fact, a RSB transition should cause a 
localization transition at the critical iFamp m- However, 
it is unclear whether or not an RSB transition would leave 
a trace in infinite-temperature ISs. 

Anharmonicity We go beyond the harmonic approx¬ 
imation, and take in account the relationship between 
different ISs. 

We study the reaction of the system to a force along a 
direction |7r), normalized to one: = 1- We exam¬ 

ine the softest mode, that is localized, and we compare it 
with the behavior of the eigenvectors in bulk of the p(A), 
that are delocalized. Therefore we choose [tt) = [tto) 
(softest mode) and | 7 r) = Ittrand), a vector whose com¬ 
ponents are chosen at random. The vector Ittrand) is 
not an eigenvector of A4, but it is a random linear com¬ 
bination of all the eigenvectors of the system. Since the 
bulk eigenvectors overwhelm the soft modes by number, 
Krand) will be representative of the bulk behavior. 

With the application of a forcing along | 7 r), Hamilto¬ 
nian 0 is modified in 

N 

^ ^ Jxy^x ' Sy ^ ^ {^x A'pTTx^ * Sx , 
||x-yi| = l a; 

( 3 ) 

where Ap is the amplitude of the forcing. If A-p > 0 
{Ap < 0), spins tilt toward (against) | 7 r). We can mea¬ 
sure quantitatively this response of the system to the 
forcing through to = J^x ^x ' ^x- We are interested in 
forcings Ap both in the linear response regime, and just 
out of it. 

We stimulate the system with forcings of increasing 
amplitude, and study when this kicks the system out of 
the original inherent structure. Ideally, the forcing am¬ 
plitude Ap would grow continuously. We simplify the 
analysis by choosing Ap = Aih, where A is a carefully 
tuned amplitude (see below) while ih is an integer. The 
unperturbed Hamiltonian corresponds to ih = 0, while 
ih = FNp is our maximum forcing (note that Fih forc¬ 
ings are not equivalent due to anharmonicities). 

This is how we check if new states were encountered 
upon increasing the forcing: (i) For each ih, start from 
the IS |s of the unperturbed Hamiltonian "Hrp. {ii) 
From minimize the energy using 'Hp{ih), and find 

a new IS for the perturbed system, |IS(i;i)). (m) From 
|IS(z?i)) minimize the energy again, using 'Hp(O) = Hrp, 
and find the IS |IS*) (with elements (iv) If 


Inherent structures that were obtained by relaxing an infinite- 
temperature configuration. 
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Figure 3. Overlap between initial and forced IS (if another 
IS is reached) as a function of the forcing step in. for H^uip = 
0.1, for forcing along the softest mode. The inset shows the 
distribution of the overlaps of randomly found ISs. 


rearrangements (at i?amp = 10, 50 the energy landscape 
is too trivial and the forcings never lead to a new IS), as 
it can be seen from Fig. top, where we show only the 
overlap with the largest forcings, of ih = 10. We plot 
1 — Qif and put it on log-scale so it is better visible. The 



I IS*) = the second minimization lead the system 

back to its original configuration, so the forcing was too 
weak to break through an energy barrier. On the con¬ 
trary, if I IS*) 7 ^ 15*0^)) the forcing was large enough for 
a hop to another valley. 

To ensure well-defined forcings along Ittrand), we nor¬ 
malize Ap with IIIt?)II \^x\- Indeed, the perturba¬ 
tion in ^ is bounded by |X)a; • 4;| < Ex ' ^xl < 

|l|7r)||j^. So, the perturbation is made extensive by choos¬ 
ing yip = Aih, with A = , where the amplitudes 

A are an external parameter (of order 1), that we tuned 
in order to be in the linear response regime for small i/i, 
and just out of it for ih approaching Np (see appendix). 

For the softest mode we analyzed the effect of forcings 
of order 0(1) because larger ones lead the system out 
of the linear response regime, so the amplitude of the 
forcings along Itto) is Ap{ih) = |yf^- 

See the appendices for further details about the linear- 
response regime, hops between valleys and the phe¬ 
nomenology of these rearrangements. 

Two-level systems In the spectrum of AI, p(A), there 
is an extensive number of very soft modes, with a lo¬ 
calized eigenstate. The eigenstates can connect different 
ISs through the forcing procedure we described. The 
connection caused by such states is privileged, because 
the couples of ISs are very similar one to the other. We 
show this in hgure where we compare the overlap qa 
between the configurations obtained through a forcing of 
amplitude Ap{ifi) with the typical overlap between inde¬ 
pendent ISs. This happens for every field that produces 


® The overlap between and |IS(*/i)) is defined as gjf = 

-k AOa, 'jif.x. with = 4)^^ • Sx{ih), being Sx{ih) the spins of 
the configuration |IS(i(,)). 


Figured. Top: 1 — gif for rearrangements that occurred at the 
10**' forcing step, for fields 4/amp = 0.05, 0.1, 0.5,1, 5. Bottom: 
cumulant W as a function of 1/N for Ittrand) and lifo). In 
both plots, the straight line is a reference curve oc 1/N. 


overlaps qi{ are much closer to 1 than the overlaps of in¬ 
dependent ISs (Fig.|^ inset). This means that the ISs are 
somewhat clustered in tiny groups that are represented 
by a single IS. This could be an operational definition of 
classical TLS, i.e. a system in which there are two very 
close states, where the transitions from one state to the 
other can be treated as independent of the rest of the 
system 

Moreover, the energy barriers separating these privi¬ 
leged states are positive, but they do not grow with the 
system size (see figure 14 in the appendix). This sug¬ 
gests that in the thermodynamic limit |IS*) and |s^^®4 
are separated by an infinitely small energy barrier, just 
as in a TLS. 

We can get more insight on the type of rearrangement 
that took place during the valley change, by defining the 


cumulant W = 


where Wx = 1 — 9if,x- If the 


rearrangement is completely localized, W = 1, whereas 
if it is maximally delocalized W = 1/N. Fig bot¬ 
tom, shows that, as expectable, the rearrangements are 
localized when we stimulate the system along the softest 
mode, and delocalized when it is along a random direc¬ 
tion. 


Conclusions The introduction of a random field in 
the Heisenberg spin glass model, besides extinguishing 
the rotational symmetry, changes qualitatively the shape 
of its DOS. Very strong random fields suppress the soft 
modes, and a gap appears in the DOS g{u}). Still, soft 
modes do resist the application of a random field when 
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it is not too large. The data are compatible with the 
absence of a gap, where for small w the DOS grows as 
g(uj) oc differently from the zero-field expectation, 
g{oj) oc [131 mi¬ 
lt appears that a finite fraction of the modes is local¬ 
ized, suggesting a localization transition when the system 
becomes large. 

We also analyzed the anharmonicity of the energy land¬ 
scape, by imposing an external force on the system. The 
reaction of the spin glass has a strong dependency on the 
direction of application of the force. Forcings along a 
random direction, need to be extensively strong in order 
to move the orientation of the spins. Equivalent results, 
instead, can be obtained through forcings of order 1, if 
they are oriented along the softest mode. 

Forcings along the softest mode cause localized rear¬ 
rangements that lead the system to a new IS that is in¬ 
finitely similar to the original one. The two states are 
separated by very small energy barriers. This could be 
used as an operational definition of classical TLSs. 
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assistance provided by BIFI-ZCAM (Universidad de 
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Appendix 

Parameters of the simulations 

In Tab. |l] we resume how many samples we simulated 
for each couple {L,H) both for the spectrum p(A), and 
for the forcings, the number of eigenvalues n\ that was 
computed with Arpack, and the amplitudes A of the forc¬ 
ings. 


Reaching the inherent structure 

As energy minimization algorithm we use the succes¬ 
sive overrelaxation, that was successfully used in [37] for 
3d Heisenberg spin glasses. It consists in an interpola¬ 
tion, through a parameter A, between a direct quench, 
that aligns all the spins to the local field hx (see. e.g. [48] 
for a recent application in Heisenberg spin glasses), and 
the microcanonical overrelaxation update (well explained 
in [33]). 


-^^amp 

L 

^samples 

^replicas 

n\ 

A(|f?o)) A(| 

ttrand)) 

50 

192 

10 (0) 

2 

35 

- 

- 

50 

96 

10 (10) 

2 

80 

- 

1 

50 

48 

70 (70) 

2 

500 

1 

1 

50 

24 

100 (100) 

2 

500 

1 

1 

50 

12 

100 (100) 

2 

500 

1 

1 

10 

192 

10 (0) 

2 

35 

- 

- 

10 

96 

10 (10) 

2 

80 

- 

0.72 

10 

48 

70 (70) 

2 

500 

0.6 

0.72 

10 

24 

100 (100) 

2 

500 

0.3 

0.72 

10 

12 

100 (100) 

2 

500 

0.3 

0.72 

5 

192 

10 (0) 

2 

35 

- 

- 

5 

96 

10 (10) 

2 

80 

- 

0.3 

5 

48 

70 (70) 

2 

500 

0.014 

0.3 

5 

24 

100 (100) 

2 

500 

0.02 

0.3 

5 

12 

100 (100) 

2 

500 

0.024 

0.3 

1 

192 

10 (0) 

2 

35 

- 

- 

1 

96 

10 (10) 

2 

80 

- 

0.05 

1 

48 

70 (70) 

2 

500 

0.004 

0.05 

1 

24 

100 (100) 

2 

500 

0.0045 

0.05 

1 

12 

100 (100) 

2 

500 

0.0045 

0.05 

0.5 

192 

10 (0) 

2 

35 

- 

- 

0.5 

96 

10 (10) 

2 

80 

- 

0.022 

0.5 

48 

70 (70) 

2 

500 

0.008 

0.02 

0.5 

24 

100 (100) 

2 

500 

0.009 

0.022 

0.5 

12 

100 (100) 

2 

500 

0.009 

0.022 

0.1 

192 

10 (0) 

2 

35 

- 

- 

0.1 

96 

10 (10) 

2 

80 

- 

0.012 

0.1 

48 

100 (70) 

2 

500 

0.006 

0.012 

0.1 

24 

100 (100) 

2 

500 

0.1 

0.012 

0.1 

12 

100 (100) 

2 

500 

0.1 

0.012 

0.05 

192 

10 (0) 

2 

35 

- 

- 

0.05 

96 

10 (10) 

2 

80 

- 

0.011 

0.05 

48 

100 (70) 

2 

500 

0.06 

0.011 

0.05 

24 

100 (100) 

2 

500 

0.42 

0.011 

0.05 

12 

100 (100) 

2 

500 

0.36 

0.011 

0.01 

192 

10 (0) 

2 

25 

- 

- 

0.01 

96 

10 (10) 

2 

80 

- 

0.016 

0.01 

48 

100 (70) 

2 

500 

0.045 

0.016 

0.01 

24 

100 (100) 

2 

500 

0.009 

0.004 

0.01 

12 

100 (100) 

2 

500 

0.007 

0.001 


Table I. Samples and replicas of our simulations. The number 
between parenthesis is the amount of samples used for the 
forcings. We indicate with n\ the number of eigenvalues we 
calculated from the bottom of the spectrum p(A). AdyfRANo)) 
and A(|7fo)) are the forcings’ parameters. 


We propose sequential single-flip updates with the rule 
^SOR ^ fex + As 

Ildx + AsOR 


(4) 
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where s is the overrelaxation update 


?OR 



7*old 


hi 


?old 


( 5 ) 


The limit A = 0 corresponds to a direct quench that 
notoriously presents convergence problems. On the other 
side, with A = oo the energy does not decrease. 

It is shown in m that the optimal value of A in terms 
of convergence speed is A « 300. Thus, the search of ISs 
was done with A = 300, under the reasonable assump¬ 
tion, that we will reinforce right away, that a change on 
A does not imply sensible changes in the observables we 
examine. In fact, the concept of IS is strictly related to 
the protocol one chooses to relax the system, and on the 
starting configuration. From m our intuition is that de¬ 
spite the ISs’ energies do depend on these two elements, 
this dependency is small and we can neglect it (depen¬ 
dencies on the correlation lengths will be examined in a 
future work [33]). 

To validate the generality of our results we compared 
the ISs reached with A = 300,1,0, at i?amp = 0 over 
a wide range of temperatures. We took advantage, for 
this comparison, of the L = 48 configurations that were 
thermalized in [^, that go from Tgc to §7 sg. 

In figure we plot the energy Eis of the reached ISs,[^ 
as a function of the temperature T. We show ten random 
samples, each minimized with A = 300,1,0. Increasing 
A the energy of the inherent structures decreases but this 
variation is smaller than the dispersion between different 
samples. The energy of the ISs also decreases with T, but 
this decrease too is smaller than the fluctuation between 
samples. Since the dispersion on the energy is dominated 


A = 300-SO - 

A = 300 - si + 

A = 300 -s2 ••><• 

A = 300 - S3 * 

A = 300 - s4 — 0 - - 
A = 300 -S5 - -■ - 
A = 300 - S6 -- -o-- - 
A = 300 - S7 
A = 300-s8 
A = 300-s9 ■* 


A = 1 - SO - 

A= 1-s1 + 

A = 1 - S2 X 

A = 1 - S3 - - 

A = 1 - s4 - -Q - 

A= 1 -S5 -- -- -- 
A= 1 -S6 - e - 
A= 1 -S7 - • - 
A = 1 - S8 A 

A = 1 - S9 * 


A= O-sO 
A= 0-Sl 
A = 0 - S2 x- • 

A= 0-s3 
A= 0-s4 -D- 

A = 0 - s5 - ■ - 

A= 0-s6 
A= 0-s7 — 

A = 0 - S8 6- 

A = 0 - S9 * 





T 


Figure 5. Energy of the inherent structure as a function of 
temperature for 10 samples chosen at random, for Flamp = 0, 
L = 48. We use the same symbol for the same sample. ISs 
obtained with A = 300 are in blue, red represents A = 1 
and purple stands for A = 0. Purple and red lines almost 
overlap. Sample-to-sample fluctuations are the largest source 
of dispersion, compared with A and T. 


by the disorder, rather than by A or T, we can think of 
putting ourselves in the most convenient situation: T = 
oo, that does not require thermalization, and A = 300, 
that yields the fastest minimization. 

Also the spectrum of the dynamical matrix does not 
show relevant signs of dependency on either T of A, as 
shown in figure]^ 



-2 0 2 4 6 8 10 12 14 16 


X 



-2 0 2 4 6 8 10 12 14 16 

X 

Figure 6. Spectrum p{X) of the Hessian matrix calculated at 
the inherent structure for Hamp = 0, L = 48, with the method 
of the moments . Top; p(A) for different temperatures from 
T — 0.12 to T = oo. Bottom: comparison of the spectrum 
between A = 1 and A = 300 at T = 0.12,0.19, oo various A. 


® With a normalization factor 1/3N that bounds it to unity. 
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The Dynamical matrix 


Derivation of the expression for the Hessian 


We calculate the dynamical matrix as the Hessian ma¬ 
trix Af of Hamiltonian 

N 

^Rf(I^) — ^ ^ Jxy^x ‘ Sy ^ ^ ' Sx ; ( 6 ) 

\x-y\ = l a; 

calculated at the local minima of the energy. It is not 
straightforward to compute At because it is necessary to 
take in account the normalization of the spins s J = 1 Va:. 

To this scope we define local perturbation vectors tTx^ 
called pions in analogy with the nonlinear a model m- 
The distinguishing feature of the pions is that they are 
orthogonal to the IS, {Sx ■ ^x) = 0, and that their global 
norm is unitary, (ttItt) = 1 , where {a\b) indicates the 
scalar product between configurations, {a\b) = ’ 

bx ■ 

We use the pions to parametrize order e perturbations 
around the IS as 

Sx = 4^®^ \/l - + CTTp, , (7) 

so the position of s' is fully determined by tTx- As long 
as e is small enough to grant < 1 Va;, the normal¬ 
ization condition is naturally satisfied without the need 
to impose any external constraint. 

We now build a local reference change. For each site 
X we define a local basis B = ei^a,,e 2 ,x|, where 

ei.x, 62 ,X are any two unitary vectors, orthogonal to each 
other and to s® , and well oriented. In our simulations 
they were generated randomly. In this basis the pions can 
be rewritten as Tf^ = ( 0 , ai, 02 ), where now they explicitly 
depend only on two components, with real values ai and 
02 . We can therefore rewrite the pions as two-component 
vectors tCx = ( 01 , 02 ). At this point we integrated the 
normalization constraint with the parametrization, and 
we can obtain the 2N x 2N Hessian matrix AI, that 
acts on 27V-component vectors Iff), by a second-order de¬ 
velopment of "His (the derivation of AI is shown in the 
subsection ahead). The obtain matrix is sparse, with 13 
non-zero elements per line (1 diagonal element, and 6 
two-component vectors for the nearest-neighbors). The 
matrix element AI“^ is 

Mll = Mxy{ea,x-ep,y), (8) 

with 

D 

f/JCis) . ?(is)\ J 5 ('91 

■‘''^xy — ’-’xyV'y ^y j / , ’Jxy'-’x+ii,y ) WJ 

fi— — D 


We derive the expression of the Hessian matrix AI of 
the Hamiltonian 71 rf [eq- (j^in the main article] that we 
implemented in our programs. 

In terms of pionic perturbations [recall eq. 0], M 
would be defined as AI“^ = ^ . An easy way 

to extract the Hessian is to write TIrf as perturbations 
around the IS and to pick only the second-order terms. 

To rewrite "Hrf as a function of the pionic perturba¬ 
tions, it is simpler to compute separately the dot prod¬ 
ucts {sx ■ Sy) and hx • Sx- Including the e factors into the 
perturbation tt^, the generic spin near the IS is expressed 
as Sx = Sx^'^ — Tr’^ + TTx- We can make a second-order 

expansion of the non-diagonal part of the Hamiltonian 
by taking the hrst-order expansion of the square root 
\/l-'^x - 1 - "^x/S, 

(Sx • Sy) = ( 10 ) 

= ( 4 ^®^ • ( 4 ^®^ + ^v) = 


= {sf^ ■ sT^) + ( 4 ^^^ • ^y) + (sT^ ■ -x) + 

+ \ (-4 - 4) (4^®^ • 4^®^) + 27fx • TTy -ko(7f3). 


On the other hand the random-held term is 
(hx ■ SxJ = hx ■ (4^®^ \/l - '^x + 


( 11 ) 


= (4 ■ 4^®^) + [hx • ff^) - Y ■ 4^®^) + o(ff^). 

By inserting eqs. ( jiTO and taking only the second-order 
terms we obtain how the Hessian matrix acts on the helds 
Itt): 


- (tT,^! AI |7fy) = 


( 12 ) 


<x,y> 

N ^2 

+27?^ ■ 4] + H Y ^ 

X 

1 ^ 

= 2 H 4^®^ • (^X ®^ +4) 




xy"y ) 


V.\x-y\ = l 


where the bold arab characters indicate the site, and 
the greek characters indicate the component of the two- 
dimensional vector. 


”* (IS) 

where we called hx the local held of the IS. The just- 
obtained expression represents a sparse matrix with a 
matrix element M.xy that comfortably splits as M.xy = 















T^xy + J^xy into a diagonal term V^y and a nearest- 
neighbor one Afxy , with 


T^xy — ^xy 


i'x , (13) 


A/*a 


xy 


= -E 

fi——d 


Jxy^i 


x+e^,y ) 


(14) 


where is the unit vector towards one of the 2d neigh¬ 
bors. 

A4 in the local reference frame The last step is to get 
an expression of the Hessian matrix in the local reference 
frame, that includes the spin normalization constraint. 

In the local reference frame the pions are written like 
if = aiei,x+a 2 e 2 ,x because they are perpendicular to the 
first vector of the basis, Sx , and that is why we write 
them in a two-dimensional representation as tt = (oi, 02 ). 

In this local basis, the matrix element acting on the 
pions is written as 


r^xd^xy'^y — '^xd^xy'^y 


(15) 


1 ■Adxy 


■Adxy (^2,x 


[ “i.y j 


' ^2,^) 

■M^xy (^2,0: 

■ ^2,y)J 

V “2,y J 


— (^1,03) ^2,ic) I 

SO in the 2A^-dimensional reference Ai is expressed as 


■^xy ~ ■Mxy {Sa,x ' ^0,y) > (16) 

and the elements of the diagonal and nearest-neighbor 
operators V and N become 


75^/3 _ X 
‘-^xy ~ ^xy 



(17) 


d 

■^xy ~ ~ ^ ^ JxySx+e^,y {^a,x ‘ ^l3,y) ■ (18) 

y——d 


Cumulatives for all the fields 






Figure 7. Cumulative distributions F[\) for small random 
fields Hamp = 0.01,0.05,0.1,0.5. In each plot we show black 
a reference curve representing the power law that is our 
guess for a universal behavior, and a grey line indicating the 
Debye behavior A^'®. One could expect a Debye behavior for 
A > A*, with A* —>■ 0 as i7amp —>■ 0. Instead, we see an 
excess of eigenvalues even compared to the Debye behavior, 
indicating a likely boson peak. Further discussions in the 
main text. 





^amp - 50 j 

L= 12 '—'— ' j 


£=24 f 


£=48 . ) 


£=96 ii 


£ = 192 



To calculate F{X), for numerical purposes we compute 
the average of the eigenvalue, and we plot k/{2N) 
versus (Afc). In this way we deliberately avoid to looking 
to the tail in F(X) that in finite volume systems is present 
as an effect of the fluctuation of the lowest eigenvalues. 
With these procedure the errorbars are on the x axis. An 
advantage is that the function F{X) does not depend on 
the number of eigenvalue computed. 

In figures and we show the function F{X) for all 
the fields we simulated. We were able to calculate with 
Arpack the lowest eigenvalues of the spectrum. The num¬ 
ber of calculated eigenvalues n\ is in table|^ All the plots 
are compared with the Debye behavior A^ ® and with the 
power law behavior A^'®, because our data suggest that 
a universality on the exponents would set them around 
around 7 = 2.5, (5 = 4 and a = 1.5. This is straight¬ 
forward for i7amp = 0.1,0.5,1,5, where when A is small 


Figure 8. Cumulative distributions F{X) for large random 
fields Ftnuip = 1,5,10, 50.In each plot we show a reference 
curve representing the power law A^'®. The orange line in the 
bottom left set is proportional to A®. 


there is a clear power law behavior, with a power close 
to 2.5, while it can be excluded for Hamp = 50, where 
the soft modes are suppressed in favor of a gap, as it was 
also visible from the inset of figure 1 of the main article. 
At Hamp = 10 we are probably close to where the gap 
forms. The F{X) goes as a large power law A® when A 
is large, but at the smallest values of A, recovered from 
L = 192, there is a slight change of power law towards 
something that could become 2.5. One could also argue 
that a F{X) goes to zero as a power law for any finite 
Hamp, as long as one looks at small enough A. Numerical 
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analysis cannot answer to questions of this type, but still, 
even if no sharp transition is present, an empirical gap is 
clearly present for large 7?amp, since the precision of any 
experiment (numerical or real) is finite. In the case of the 
smallest fields i?amp = 0.01,0.05, we suffer from effects 
from i?amp = 0. The spin waves do not hybridize with 
the bulk of the spectrum, and pseudo-Goldstone modes 
with a very small eigenvalue appear, making it hard to 
extract a power law behavior. 

Overall, we see good evidence for a 7 around 2.5 at 
several values of I?amp, and at other fields the data is not 
in contradiction with a hypothesis of universality in the 
exponents defined by 

p(A) ~ A“ , F{X) ^ . (19) 


Making the reasonable assumption that different eigen¬ 
vectors do not interfere with each other, and exploiting 
the orthogonality condition '^^ipm(x)ipn(x) = 5mm we 
obtain the desired correlation function 

C{r) =g^{r) = + . (23) 

This correlation function favors the softest modes by a 
factor 1/A^. This is an advantage, because the bulk 
modes do not exhibit a finite correlation length, so it 
is useful to have them suppressed. 

The correlation function for small fields 


When the field is small we remark a change of trend from 
7 « 2.5 to 7 < 1.5 at a value A*. The crossover A* shifts 
towards zero as i?amp decreases. This probably indicates 
the presence of a boson peak, an excess of modes at low 
frequency. Signs of a boson peak in at i?amp = 0 can be 
seen in figure In that case the mass of the spectrum is 
all concentrated at low A, but there ought to be a Debye 
behavior, meaning that A* is very little. 


Eigenvector correlation function 

Since in a localized state the eigenvectors have a well- 
defined correlation length, we can use also this criterion 
to probe the localization. We can define a correlation 
length from Green’s function Q, that is defined through 
the relation JHQ = S^y, an is commonly used in field the¬ 
ory for two-point correlations. Since shares eigen¬ 

vectors |'0„) with Ai and has inverse eigenvalues 1/A„, 
Green’s function is [3 

( 20 ) 

n 


and squaring the relation 


As we show in Fig. when i?amp is small the correla¬ 
tion function is not exponential. 



r 


Figure 9. The correlation function C(r) extracted from the 
eigenvectors in L = 192 lattices. 


Forcings 




XmXn 


( 21 ) 


By averaging over the disorder we gain translational in¬ 
variance and can be written as a function of the dis¬ 
tance r = X — y, 




[‘)pm{x)tpn{x)][lpm{x + r)lpn{x + r)] 

V ■ 


( 22 ) 


^ For simplicity we use Ai-component eigenvectors ipnix) instead 
of the 2Ai-component ones lir). The relationship between the two 
can be recovered through 1 / 1 ^ (cc) = 


Probing the linear regime 


To make sure that our forcings are not too strong, we 
monitor the direct reaction of the system to the forcing. 
We define a “polarized magnetization” m = (IS(z?i) | 7 r) = 
that indicates how much the forcing pushed 
the alignment of the spins along the pion. The amplitude 
of the forcing is tuned well if m(*;i) is close to the linear 
regime. In table |T] we show the amplitudes A we used in 
order to be in the linear regime. Figure 10 confirms that 
this was the working condition for the forcings along [tto). 
Figure 11 is analogous, but along Ittrand)- In the latter 
figure we rescale to by a factor l/^N to obtain a col¬ 
lapse. In fact the normalization (ttrand Ittrand) = 1 
implies that the components of Ittrand) are of order 
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1/Vn, so the polarized magnetization is bounded by 
|m| = \{IS{ih) Ittrand) I < See ~ 


4 



0.5 

^amp = ^ T ’ 



0.25 


S 0 



<S 

0 


-4 


■■■ . . . 

-0.25 

Z,= 12 - 

• ■ . ‘ Z, = 24 



•0.0005 0.0005 

-0.5 

.i= 48 _ . 


■0.03 -0.02 -0.01 0 0.01 0.02 0.03 -0.002 0 0.002 




Since the spins in our model are continuous variables, we 
impose AW = 1/lV as a threshold to state whether there 
was or not a change of valley. 

The cumulant W is an indicator of the type of rear¬ 
rangement that took place. If the rearrangement is com¬ 
pletely localized (only one spin changes), W = 1, whereas 
if it is maximally delocalized (all the spins have the same 
variation), then W = 1/A^. 


Figure 10. Polarized magnetization m of the forcings along 
Itto), for Jfamp = 0.1 (left) and Jfamp = 1 (right). The inset 
is a zoom of the same data. 




‘h 


•h 


Figure 11. Rescaled polarized magnetization m of the forcings 
along Ittrand), for iFamp = 0.1 (left) and iFamp = 1 (right). 
The data are rescaled in order to collapse. 


Ending in a new valley. For each Ap{ih) we measure 
the overlap qif between the two minimas of Hrf, the ini¬ 
tial IS, and the final one, |IS*). Naively, checking 

that gif < 1 in principle is a good criterion to establish 
whether the system escaped to another valley. We pro¬ 
ceeded similarly, in terms of the spin variations between 
initial and final configuration, through the quantities 


N 

|is*) = iV(l - gif), 

X 



(24) 

(25) 

(26) 


The local variation Wx measures the change between the 
beginning and the end of the process. If the spin stayed 
the same then Wx = 0, while if it became uncorrelated 
with the initial position Wx = 1 in average. If one and 
only one spin becomes uncorrelated with its initial config¬ 
uration, the variation of W is AW — 1/N. Similar vari¬ 
ations AW do not mean that one spin has decorrelated 
and the others have stayed the same, this is impossible 
because and |IS*) are ISs and collective rearrange¬ 

ments are needed. A AW = 1/N means instead that 
the overall change is equivalent to a single spin becoming 
independent of its initial state. This is, for a rearrange¬ 
ment, the minimal change in the W that we can define. 


Falling back in the same valley. Even though the 
forcing is along a definite direction, since the energy land¬ 
scape is very irregular, it may happen that stronger forc¬ 
ings lead the system to the original valley. For example 
it may happen that ih = 2 lead the system to a new val¬ 
ley, and ifi = i lead it once again to the same valley of 
ih = 1. To exclude these extra apparent valleys we la¬ 
bel each visited valley with its W, and assume that two 
valleys with the same label are the same valley. These 
events are not probable, and even less likely it is that this 
happen with two different but equally-labelled valleys, so 
we neglect the bias due to this unlucky possibility. 


Rearrangements 


To delineate the shape of the energy landscape, we 
want to study, for every couple (i?amp,T), the probabil¬ 
ity that a forcing of amplitude Ap lead the system to a 
new valley, to distinguish the behavior of soft from bulk 
modes. 

Furthermore, once the system made its first jump to a 
new valley, it is not excluded that a bigger forcing lead it 
to a third minimum of the energy. One can ask himself 
what is the probability Ph^^^,l{Ap^ n) that n new valleys 
are reached by forcing the system with an amplitude up 
to Ap{ih), and to try to evince a dependency on sistem 
size and random field. Even though n is bounded by z^,, 
this does not necessarily mean that if we made smaller 
and more numerous forcings n could not be larger. On 
another side, if for a certain parameter choice rearrange¬ 
ments are measured only for large ih, it is reasonable to 
think that these represent the smallest possible forcings 
to fall off the IS. 

To construct PH^^j,,LiAp,n), for every replica and 
sample we start from ih = 0 and increase \ih\ either in the 
positive or negative direction (the two are accounted for 
independently). The value we assign to PH^^p,L{Ap,n) 
is the number of systems that had n rearrangements after 
ih steps, divided by the total number of forcings, that is 
2N N 

’ rep-^ ’sam • 

First rearrangement. In figure we show the prob¬ 
ability of measuring exactly n rearrangements after ih = 
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TVp = 10 forcing steps. Even thongh both for Ittrand) 
and Itto) we are in the linear response regime, the behav¬ 
ior is very different between the two types of forcing. In 
the first case every single forcing step we impose leads the 
system to a new valley. In the second rearrangements are 
so uncommon that even though the probability of hav¬ 
ing exactly one rearrangement is finite, that of having 
more than one becomes negligible for large samples. It 




Figure 12. Probability of there being exactly n changes 
of valley after ih = A^p = 10 forcing steps. The data 
come from //amp = 0.1 (left) and = 1 (right). If 

PH^^p,L{AF,n) = 1 for n = 0 it means that the forcings were 
not strong enough to ever get out of the initial IS. On the 
contrary, /’namp,n(^F, n) = 1 for n = 10 means that every 
single step lead the system to a new IS. The latter scenario is 
realized in the case of forcings along Ittrand), especially when 
the system size is large. On the other side, forcings along Itto) 
display a small but finite amount of rearrangements. 

is then reasonable to think that any rearrangement we 
measure for Itto), it occurs for the smallest possible forc¬ 
ing, and even when more than one occurs, these jumps 
are between neighboring valleys, where by neighboring we 
mean that no smaller forcing would lead the system to 
a different IS. To convince ourselves of this we can give 
a look at the average number of rearrangements after ih 
forcing steps, n{ih) (figure [l^ . When ih is small no 
new ISs are visited and n(ih) = 0, while for larger ih, 
n{ih) is positive but small, so we can call these changes 
of valley “first rearrangements”, i.e. rearrangement be¬ 
tween neighboring valleys. 


Energy Barriers 

To reinforce the idea of two-level system (TLS), we 
report the energy barriers between the couples of con¬ 
nected ISs. The maximum value of AE before a hop 
to another valley should give an estimate of height of 


® We do not show data regarding forcings for i/amp = 10, 50, be¬ 
cause no arrangement takes place. Most likely the energy land¬ 
scape is too trivial. 

® Because (Ap, n) is not defined over all the samples (it is 

hard to reach many different valleys and it may not happen in all 
the simulations), the errors on (Ap,n) were calculated 

by resampling over the reduced data sets with the bootstrap 
method. 




'ft ‘A 


Figure 13. Average number of rearrangements n{ih) for forc¬ 
ings along Itto) (left) and along Ittrand) (right). The data 
come from //amp = 0.1 (top) and //amp = 1 (bottom). When 
the lattice becomes large enough, the forcings along Ittrand) 
lead to a new IS every time ih is increased. The data from 
the |ffo) and //amp = 1 can be said to be in the regime of first 
rearrangement. 


the barrier. Still, it may happen that the IS obtained 
with the Hamiltonian T-Lp have an energy lower than the 
energy Erf( |s^^®^) ) calculated with "Hrf, so in a strict 
sense AE is not positive definite. To overcome this issue, 
we measure the barrier AE* from the arriving IS instead 
of the starting one, that has the advantage of being pos¬ 
itive definite. We see in figure [T^ that while the energy 
barriers from random forcings increase with the system 
size (the growth is 0{N)), while those within the TLS 
(along the softest mode) do not grow at all. 



Figure 14. Average energy barrier AE* for forcings along 
Ittrand) (top) and |ffo) (bottom), for random fields of ampli¬ 
tude //amp = 0.1 (data for |ffo) and L = 96 were not com¬ 
puted, see table I). 
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